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CONCEPTS FOR RADICALLY INCREASING THE NUMERICAL 
CONVERGENCE RATE OF THE EULER EQUATIONS 


I . ) Introduction 

There is considerable interest at the present time in 
solving the Euler equations to predict fluid flow. The Euler 
equations are assumed to be good approximations of the Reynolds 
averaged Navier-Stokes equations when modeling attached flows 
with shock waves and rotational effects as well as being much 
cheaper to compute because the fine details of the physical 
viscous region are neglected. 

All of the current methods of solving the Euler equations, 
for example those of Reference 1, rely on iteration procedures 
that linearize the governing equations about values at the 
previous iteration level. The impetus behind the present study 
is to determine a better function or base solution for the 
linearization in order to reduce the number of iterations and, in 
particular, to determine whether the number of iterations can be 
reduced to one. 


In Reference 2, which is the report on the Phase I effort, a 
study of the transonic small disturbance (TSD) equation is 
described. The idea behind the study of the TSD equation is that 
it is a simple nonlinear equation which can represent shock 
waves. If this equation can be linearized so that shock waves 
can be represented then this technique may be transferable to the 
Euler equations. In Reference 2 applications of the transonic 
perturbation technique, developed by Nixon (Ref. 3), are 
investigated with the conclusion that a linearized solution of 
the TSD equation is feasible. However, during the course of the 
study it was found that a solution of the linearized partial 
differential equation obtained by the integral equation technique 
differed greatly from solutions obtained by a finite difference 
technique. Since it is anticipated that general solutions to the 
linearized Euler and TSD equations will be obtained by finite 
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difference methods, this is a disturbing result. It is suggested 
also in Reference 2 that a closed form solution of the TSD 
equation could be found. 

This report is concerned with extensions of the studies 
given in Reference 2, in particular a more detailed examination 
of the TSD equation and the extension to the Euler equations. It 
had been the intention to extend the investigation to include 
three dimensional flows but lack of time and funds precluded this 
step. The report is written as a series of relatively self 
contained sections to facilitate understanding. In Section (2) 
the anomaly between integral equation and finite difference 
solutions of the TSD equation reported in Reference 2 is 
resolved. The conclusion is that the finite difference technique 
applied a fictitious constraint on the solution thus leading to 
erroneous results. 

Section (3) is concerned mainly with further studies of the 
integral equation formulation of the TSD equation. One of the 
more surprising results of the study is that the original 
formulation of the transonic perturbation method given in 
References 3 and 4 is incorrect although the published results 
are correct because of a numerical smoothing in the computations . 
An important result is that an arbitrary solution of the integral 
form of the TSD equation, that is a solution that does not locate 
the shock wave, gives oscillatory behavior in the neighborhood of 
the sonic line. This oscillatory behavior is a numerical 
representation of the mathematically derived expansion shock. 
The correct shock location is that which removes these 
oscillations, which indicates that the oscillations are a 
critical consequence of the formulation. A dissipative finite 
difference solution of the linearized TSD equation removes these 
oscillations by numerical means and thus removes the necessary 
information for determining the shock location. A non- 
dissipative finite difference solution is derived which does 
locate the shock wave correctly. The results of both integral 
equation and finite difference results indicate that 
linearization of the TSD equation can lead to adequate results . 
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In Section (4) an attempt to use a linearized TSD equation 
with constant coefficients is described. If sufficiently 
accurate, this would allow a transonic solution to be constructed 
using superposition of two Prandtl Glauert solutions, one 
subsonic and one supersonic. The results of this study are 
inconclusive; good results can be obtained but the technique 
seems to be very sensitive to numerical errors . 

In Section (5) the techniques developed for the finite 
difference solution of the TSD equation are applied to the 
linearized Euler equations. It is found that the odd and even 
points in the Euler solution have to be decoupled before adequate 
solutions can be obtained. It is concluded that solutions of the 
linearized Euler equations can be of adequate accuracy, although 
it must be stated that this conclusion is based on limited 
experience . 

In Section (6) a study is made of separation phenomena 
modeled by the Euler equations . It is concluded that although 
separation can be modeled by the Euler equations, it does not 
have much in common with real (physical) viscosity controlled 
separation. It is found also that specifying an empirically 
determined separation line is not consistent with the Euler 
equation formulation. Finally, it is suggested the boundary 
conditions commonly used in numerical solutions of the Euler 
equations are inconsistent. 


2 . ) A Study of a Discrepancy Noted in Phase I 

In the Phase I work two methods of calculating the 
linearized TSD equation were used, namely an integral equation 
method and a finite difference method. The TSD equation is 
linearized about a piecewise constant base flow with one portion 
giving an elliptic equation and the other a hyperbolic equation. 
In the finite difference solution, calculated using a modified 
version of the TSFOIL code (Ref. 5), the solution in the 
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hyberbolic zone is subsonic, while part of the flow in the 
elliptic domain is supersonic. This anomalous behavior must be 
understood before further work is performed since it may be due 
to inherent difficulties in the finite difference approximation. 
An example of the anomalous behavior is shown in Figure 1. 

The computational code used in these computations is TSFOIL 
(Ref. 5) which solves the equation 

(1-u )<b + <p = 0 (1) 

v o /Y xx Y yy v ' 


where <p is a perturbation velocity potential and u Q = 0 or 2; 
when u Q = 0 the equation is elliptic and when u Q = 2 the equation 
is hyperbolic. At the "sonic line" (upstream limit of the 
hyperbolic domain) the algorithm in TSFOIL puts 



( 2 ) 


since in a continuous distribution of u Q , u q is unity (sonic) at 
this boundary. However, since u Q is discontinuous in Equation 
(1), the application of Equation (2) makes 


<t> 


xx 


8u 

8x 


0 


(3) 


at the upstream limit of the hyperbolic zone. 

The differential equation will admit discontinuous solutions 
of the form 


+ <P X = 0 (4) 

for a discontinuity normal to the x axis. Since the TSFOIL 
algorithm forces <p to be zero at the sonic line it follows that 

XX 


- 6 - 


( 5 ) 



Thus the only compatible solution is 

*x = tx = 0 ( 6 > 

and the flow in at least the initial part of the hyperbolic 
domain is subsonic ( <p • < 1). If the flow accelerates over the 

succeeding portion of the airfoil, as is usual, then it is 
possible that at the downstream limit of the hyperbolic domain 
the velocity is still subsonic, in which case an expansion shock 
wave must exist there. This is the situation in the computations 
shown in Figure 1 . The anomaly is due to the enforcement of a 
continuous <p x at the "sonic" line, a condition which is not a 
solution of the differential equation. 


3 . ) Studies of the TSD Equation 

3.1 Introduction 

The TSD equation contains several of the features of the 
Euler equations, in particular the presence of shock waves. The 
equation is relatively easy to analyze compared to the Euler 
equations and will be used to provide insight into some of the 
problems arising in the linearization. 

The mechanics of correctly linearizing the TSD equation are 
examined in detail in Reference 3. Basically the nonlinear TSD 
equation is expanded in a power series in a small parameter, e, 
in terms of a coordinate system that is also perturbed as a 
series. This strained coordinate technique allows a valid 
perturbation even when shock waves move during the perturbation. 
The original work was for steady flow; a later extension is to 
unsteady flow (Ref. 4). The objective of this work was to 
perturb the airfoil geometry or the flow parameters to obtain a 
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solution in the neighborhood of the original base solution. The 
objective in the present work is to allow large perturbations 
from the base solutions and thus estimate the accuracy of the 
linearized TSD equation. The base solution need not be related 
to the problem under consideration other than by the fact that it 
has a supersonic zone with a shock wave. It is pointed out in 
Reference 3 that if the base solution does not have a shock wave 
then the linearized equation cannot have a shock wave. 

The work reported in References 3 and 4 uses an integral 
equation analysis which allows a coupling of the boundary 
conditions and the velocities to form two coupled integral 
equations. These equations are easier to analyze than the basic 
partial differential equations and their associated boundary 
conditions. However, since the object of the present study is to 
extend the linearization techniques to the Euler equations, it is 
desirable to use finite difference methods because of their 
advanced state of development. Accordingly the study described 
in this section starts with a more detailed examination of the 
integral equation than that given in Reference 3 and then extends 
the basic ideas to finite difference methods. 

The first result to emerge from the study is that the 
perturbation theory given in References 3 and 4 is in error, 
although the results are correct. The reasons are given later in 
this section. It is found that the shock location is such that 
oscillations in the neighborhood of the sonic line vanish. This 
criteria, found using the integral equation, applies to the 
finite difference formulation. As a consequence it is essential 
that the finite difference approximation does not, in general, 
remove these oscillations since their presence is critical to a 
correct location of the shock wave. In this section, results 
found using mixed differencing are incorrect because the 
difference scheme removes these oscillations, thus allowing the 
shock wave to be located anywhere. 
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Both the integral equation and the finite difference 
approximations are used to compute solutions of the linearized 
TSD equation, and these results are of adequate accuracy. 

3.2 Basic Equations 

The basic transonic small disturbance equation is 

^xx + ^yy - ^x ^xx 

with the tangency boundary condition (for a symmetric airfoil) 

</>y <x, ± 0) = y^(x) (8) 

The far field boundary condition is that <j> , <p + 0 at an 

^ r y 

infinite distance from the airfoil. Change the coordinates to 
(x' ,y ' ) where 


x' = x + f(x,y) 

y' = y + g<x,y) 


( 9 ) 


The differential equation. Equation (7), can then be written 


as 


^x'x' + ^y'y' ^x'^x'x' + 2 9x' + ^x^ ^x' + g x ^y'^ ^x'^ 


+ f x a! 7 {[(1 + f x>*x' + g x ^y'l 2/ 2 " [(1 + + ^ K'U 


x r x 3 x r y ' 


+ g x ar {C(1 + f x^x' + g x ^y'3 2/ 2 - + + 


X r X 


3 x r y i 


f v 3x' + g y^ ^y' + f y^y'^ g y 9y' + g y^y' + f y^x^ 


y ax' 

9 


a 


9x' ^x^x' + g x^y' ^ 9y' ^ g y ^y' + ^y^x^ 


( 10 ) 
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and the boundary condition becomes 


4> y (*, ± 0) = (1 + g y ) <f> yl + f y ^ x , = y' s (x(x')) (11) 


Denote 


(1 + f x^x' + g x^y' by U 


(1 + g y )^y» + f y0 x ' b Y V 


Equations (10) and (11) can then be written in integral form 
using Green's identity to give 


TT u? i_ f i Aty>(C(^)) - - £) 

2 2n i , (1 2 . ...2 


(x' - C ') 2 + y ' 2 


. U 


- / S / K C'X' tj- - U**' ~ Vn' }dS 

f ' 1 2 - -J ’ *£ at)' l ~'2 " WJ " 9 rj- ^rj^rj' 


+ / s / K x- < f f at 7 if” " u) + g e an 7 t u/ 2 - - S^ 7 + tj e . 1 


- f V - g V}dS 

rj 3f' *rj dr]' 1 


’*rs 

( 12 ) 


where 


K (x, f ; y,» 7 ) = ^ In {(x - f) 2 + (Y ~ V) 2 ) 


(13) 


and S is the flow domain shown in Figure 2. Now let 


U = U + U, 
O 1 


V - V o + V 1 


(14) 
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where U Q and V Q are some specified velocities: Equation (12) then 
becomes 

Jj2 2 

“l' 1 - °o> - 2 1 ■ 7 L - ^ K J'x' (U o D l + I 1 - A > ds 

+ Wt f £ ef 7 + g i 5r lt(U o + D l> - (U o + D l> 2 / 2 ])dS 

- V K x'< f , af 7 (v o + V 1> + 9 t W (v o + V >ds 

- W<a!r c,ds <15) 

where 


I L 2tt I 


1 i A[y' (£<£')) - f B] (x - 


(x' - C') 2 + y' 2 


d^' 


+ V - u o - V K {'x' u o' 2 ds 


(16) 


and 


B 


. < D o i °!> t f f + f f g n - 2i V * g f (V ° + v i ) 

(d + g,xi * f £ ) - f,g { J 

[<1 + g,)(U 0 + Uj) - g ? (V D + V 1 )]/[(l + f £ )(l + g^) - f^J 


(17) 


C - [(l + £,)(V o 


+ V - V" O + u i>l/t<i + f { )(! + 9 n ) - *,g £ J 


The equations given above are exact; no approximations have 
been made. The velocities U , V Q are arbitrary. As a first 
approximation let 


'V < < IU o l 

lv x l < < lv o l 


(18) 
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such that the nonlinear terms in or V.^ can be neglected. 

Also, assume that products of and V 1 with derivatives of f and 
g can be neglected. This leads to the approximate form of the 
equation 

V 1 - U o> - : L - - * ,dS 

+ J s S K x'<t f ( bI 7 + 9 ( W 1 [U ° 

aV o 9V o 

’ ^ 5^' ^rj dr}' ^ dS 


- V 2 1 - 


ac 

drj’ 


( 19 ) 


and 


A 


U ff_ + f„q - q^f 1 + q^V 
o l f tL. l£ 0 

[(1 + g^) ( 1 + f^> - f^] 


(20) 


c - t v o % + f cS “ W + Vo ]/t(1 + V (1 + V " Vc ] 


In the general cases discussed later the y straining is not 
a function of x and hence 


g^ = 0 


( 21 ) 


This leads to the simplified equation 


V 1 - V - : L - VVx'fVl - A > ds 


+ JJK Atf S77 [ 0 „ - v 2 / ,] 



and 


A 


U f.(l 

° i_ 


i? 


i(i + 


V (1 


V’ 


c 


tw 1 + V + w 

[1 + g^)(i + f € >] 


(23) 


I T is the same as I T defined in Equation (16) but with B replaced 
by B. 

where 


B - U Q / (1 + f^) 


(24) 


The next step in simplification is that the coordinate 
straining terms f, g are small and expressions containing them 
can be linearized. Thus, Equation (22) reduces to 


V 1 - V 


af 7 <u ° 

av av 

f tTFT ~ 9 a — ?}9S 
V ^rj drj' 1 


f c V ds 

U o // 2^ - dr]' ^r/ V 


a 


ar<Vo ) 

(25) 


and 



A[y'(£') + £ y; (£' ) 1 ( x> 

[(x' - £') 2 + y' 2 ] 



(26) 


Equation (25) is further simplified by neglecting f^ and 
putting g zero. The first of these assumptions implies that the 
straining does not vary significantly with y; the second 
indicates there is no straining in the y direction. 
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Neglecting and g # Equation (25) can be written as 


V 1 - V ■ J L1 - + 5 Vf 


( 27 ) 


where 


l_ , Ay 4 (D(x> - £) 

L1 2 * ° t(x' - i ’) 2 + y' 2 ] 


(28) 


z i ’ ■fs/ K ? 'x' £ «- u o ds + V K x' f ?' ax 7 (u o - u o 2 / 2> ds 


t 1 Af y- (£')(x' - £') 

+ 2V lo •> 9 


,-f ' * 2 + v' 2 


(29) 




and 


f = f/5x s 


(30) 


where 5x g is some parameter related to the shock location. 

3.3 Analysis 

Equation (27) can be written in the form 

U 1 < x - U o> * 7 III + T S 1 + Sx s r £ < 31 > 

where I is the double integral in Equation (27). This is an 
s 

equation for the unknowns and 6x g . If the integrals in 

Equation (31) are discretized then, if there are N x points in the 

x direction and N y points in the y direction. Equation (31) 

becomes a set of algebraic equations for N T values of U 1 and 

N values of 6x where 
s s 
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( 32 ) 


N = N N 
T X y 

and N is the number of y grid lines covering the supersonic zone 
in U Q . It can be seen therefore that if there is a supersonic 
zone in U Q then there is an insufficient number of equations to 
solve for the N T + N g unknowns. The additional equations are 
obtained by differentiating Equation (31) with respect to x to 
give 

a 9l f 

U lx - V - U ox U l ■ 55 < J L1 + 'ali + Sx s 55 " < 33 > 


If U lx is finite when U Q = 
expansion shock, then at x 
(33) gives 


* 


1, that is does not constitute an 
(the value of x when U = 1) Equation 


* 

- u 

ox 


★ 



ax {I li + I si }l x=x* 


6 1 4 


+ Sxl 


s^ax ^ x=x* 


(34) 


If the integrals are discretized as before then Equation (34) 
gives N equations for the N values of 5x„ in terms of U- . 

s S S X 

Hence, Equations (31) and (34) are sufficient to solve for U 1 and 
Sx s . 


The above description of the shock locating mechanism is 
adequate for the integral equation formulation but does not 
really provide much insight into how the shock is located in a 
finite difference formulation of the TSD equation. Accordingly 
some further investigation into the nature of the mechanism of 
the shock location is necessary. 

Equation (31) is an equation for U 1 and, if the integrals 
are discretized, then the following matrix equation results 



U 




B. + 5x a C. 

X S X 


(35) 
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where B. is a matrix of I T . and C. the matrix of I- at the “i"th. 

X Lil X X 

grid point and is the value of U 1 at the jth grid point. 

Equation (35) can be inverted to give 


U 






c i) 


(36) 


and some criterion must be used to select 5x„ . 

s 

In the integral formulation. Equation (34) can be used in 
conjunction with Equation (31) to give a modified equation 


A ij U lj ~ B i 


(37) 


where is a representation of the term 


* * 


V 1 - V + M K {'x' U o U 1 dS - < U ! U o - ^S^ K «-x'x- D o D l dS ) >V r f , 


and B^ is a representation of the term 


B i - <1 X . V’v 


It would appear that the use of Equation (36) would give 
some insight into the nature of the shock fitting. Consequently 
the matrix was assembled and each of the terms on the right 

hand side computed. 

In order to allow for a variable shock movement the matrix 

/ k) 

C. is split into components C^ v ' where 


C = E S C^) 
c i k-1 i 


( 38 ) 
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and N g is the number of grid points- in the supersonic domain. 
The matrix contains only those values of on the kth grid 
line with the rest of its elements zero. When k is equal to N g 
the matrix s ^ contains all of the elements in the outer 
subsonic domain in addition to those on the N g th y grid line. 
Thus Equation (36) can be written as 


In Figure 
oscillations . 
the shock wave. 


'lj = 

A Ij B i 

N 

+ k£i 

* X s k A ij 

c (k) 

c i 

(39) 

the 

term 

A ij 

is 

shown 

and exhibits 


This is equivalent to a solution without moving 


If Equation (37) is solved, the result for U, . is smooth, as 

s : « J- 3 

shown in Figure 4. If the values of s^ calculated during this 
solution are used in conjunction with Equation (39) then the term 


N 


kEi 




oscillates as is shown in Figure 5 
terras on the right hand side of 
shown in Figure 4. 


However 

Equation 


, the sum of both the 
(39) leads to the 


From the preceding discussion it follows that the correct 
shock location is the location that removes the oscillations. An 
alternative way of locating the shock wave is to remove the 
oscillation in the solution of Equation (39). One method of 
achieving this is to enforce the condition 



( 40 ) 
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on each y grid line that has a sonic point, x*. This condition 
forces the solution U in the neighborhood of x* to have curvature 
of constant sign, thus disallowing the oscillatory behavior. The 
example shown in Figure 4 was recalculated using Equation (40) 
rather than Equation (34) to locate the shock wave. The results 
of both methods were almost identical. 

Some test cases were computed using Equations (31) and (34) 
as a basis. In all the cases the base solution is a 10% biconvex 
airfoil at = 0.808. A single straining function of the type 

given in Reference 3 is used. In Figure 6 the pressure 
distribution around a NACA 0010 airfoil at = 0.808 is shown, 
and it may be seen that the agreement with the direct solution is 
fair. However, considering the magnitude of the perturbation 
from the base solution, the agreement with the direct result is 
more startling. A similar remark can be made about the result 
shown in Figure 7 where the pressure distribution around a NACA 
64A006 airfoil at = 0.875 is shown. 

In the original papers on the strained coordinate method. 
Equation (34) and its unsteady flow equivalent is not used. 
However, the solution around the sonic point was obtained by a 
second order interpolation using neighboring points. This 
interpolation actually enforces Equation (40) and hence the 
results of References 3 and 4 are correct, although the theory 
given in their papers is incomplete since it does not explicitly 
discuss the use of a condition like that of Equation (34) or 
Equation (40). In the following discussion the causes of these 
oscillations are examined. 

3.4 Origin of the Oscillatory Behavior 

If there is no shock term, that is 5x g is zero, then 
Equation (31) can be discretized on a uniform grid to give 

n t 

U i <! - U oi> = j£l f ij U ij A S + B i < 41 > 
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where is the value of I Li at the grid point i and A£ is an 
influence function that gives the contribution of Uj at the jth 
point to the ith point. 

Velocity perturbations at points i and i + 1 on a uniformly 
partitioned grid are given as 


N, 


T 


“it * 1 - °oi> - f i,i A ?l ■ jSi f i'j u j A < + f i,i+l u i+i A £ + B i 

j^i/i + l 


( 42 ) 


N 


"i+lK 1 * D o.„) - “ jfil f i'j u j A « + f i+l,i °i 4 « +B i+l 

<43) 


If AU = U^ +1 - u j_/ then AU can be expressed as 


AU 


. A t A ^ f i.i.i- f i.i+i>-< c i+i- c i>i- A ^Wi.m- B i f n-i.i> +c i B i 

( c i c i + i- f i,i + i J i + i,i< A «> 2 ) 


where 


( 44 ) 


C i 


d - U o.) 

l 


f . . Af and A 

i,x 


N 

jfil J 
j^i/ i + l 




U. AC- 


Expanding all variables at the point i+1 in terms of their value 
at the point i, using a Taylor series expansion, and neglecting 
the terms with second or higher order of A£ the following 
relation is obtained; 


AU 


AU 


AB 


< A+B i> T?li A * + ( 1 -V)ac li A * 


2 AU o 

(l-u y - (l-u ) —~ 

v o^' v o^' Af 


i A * 


( 45 ) 


- 19 - 



Next consider Equation (42) and note that 


= A + tf i,i u i + f i, i+ i< u i + ff li a £>j A * + B i < 46 > 


To the lowest order, consistent with the order of Equation (45), 
Equation (46) can be reduced to 


A + B i“ UjU - ° oi > - («i,i + f i(i+1 > “i 


(47) 


Substituting this relation in Equation (45) and neglecting terms 
2 

of (A£) the following equation is obtained. 


AU = U.^., - U. 

l+l l 


, AB A* 

+ 4 ? 


U. — T 7 . A£ + T? 
i A£ i s A£ : 

AU - 

( 1 -U ) - — 7 ^ . A£ 

v o^' A£li ^ 


(48) 


A similar procedure can be used to show that 


AU = U ± - U i ;L 


u iT7 + H 

AU | 

<!-%> + -Afli A f 


(49) 


and it can be seen that AU + # -AU . If U o^ ^ 1 then AU + and AU’ 

ft TT 

approach zero as A£ •+ o and hence 1^1 is finite. 


U o. = X ' 
i 


( 50 ) 


that is, at the sonic point. 


- 20 - 


then 


AU = -AU = U. + 


j> _ i AU . 

/A® I \ / / — ® I ' 

( A£li )/( A^ I i 


( 51 ) 


Hence U oscillates about the sonic point since 


U i-1 - U i = U i+1 - U i 


Furthermore, as A( ■> o, AU 5 * o and hence 


< 22 , *. » 

v o£ ' sonic 


( 52 ) 


which is an expansion shock. 

A measure of the amplitude of the oscillation is given by 
the second derivative, which gives in discretized form 


. AU . AU . . 

au + - au- - 2 (Tfli) + zfl'i 1 Ar 


AU 


( 53 ) 


o. „ 2 




AU .. 

If — and are considered constant then the amplitude decays 


as 


(^o.) 


increases. In other words, the oscillations decrease 


the further the point is from the sonic point. 

The conclusions reached by this analysis are consistent with 
the computational evidence shown earlier. 
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3.5 Multiple Coordinate Straining 

It had been intended to use a multiple straining of the 

coordinates so that both the shock and the sonic line location 

could be changed; the sonic line location can be fixed at the 

location of the sonic line in the base solution U . This 

o 

multiple straining would ensure that the hyperbolic domain in 

both the perturbed and base solutions would be identical in the 

strained coordinate system. Preliminary tests of this idea did 

not prove fruitful, and a modified version of the idea was 

investigated. Briefly, this modification consisted of evaluating 

the second straining parameter to minimize a weighted sum of the 

|U U I over the solution. The rationale is that since the lU^I 

is a minimum in some average sense, then the error in the 

2 

linearization process, which is of the order of lU-^l, would be a 
minimum, thus improving the overall accuracy of the solution. A 
preliminary example is shown in Figure 8 where the pressure 
distribution over a NACA0010 airfoil at a Mach number of 0.808 

is shown. It can be seen that the effect of the second straining 

is considerable. However, this is a preliminary result and 
further work is necessary to consolidate the concept of using an 
extra straining to improve the accuracy of the solution. 

3 . 6 Summary 

In this section the concept of using the linearized TSD 
equation to model transonic flow is examined. The method does 
appear feasible. The shock wave is located by requiring that no 
expansion shocks can occur in the solution, and two methods of 
enforcing this condition with the integral equation method are 
derived. One form of this condition. Equation (40), can be used 
in conjunction with a finite difference solution of the TSD 
equation in differential form. 
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3.7 Linearized Solutions of the Differential TSD Equation 

In its simplest form the linearized TSD equation, equivalent 
to Equation (27), can be written as 


*lx-x- - D o /lx- + * lyy " 2 fe- ”(1 + f x )* ox .] 2 - *o X -> 


# 2 

. f 9 ,lox' ! 

x 5x' ^ 2 ^ox' ^ 


(54) 


with the tangency boundary condition 


0 ly ( x '/ *°> = Yg ( x ( x ')) - Yg 0 ( x ' ) 

= y; (x') - y^ 0 (x') - f y' (x') 


(55) 


To investigate the solution of this equation the computer 
code TSFOIL (Ref. 5) is used. TSFOIL uses mixed upwind and 
central differencing to solve the nonlinear TSD equation. 
Equation ( 1 ) , and special parabolic and shock point operators to 
allow a conservative transfer from one scheme to the other. The 
parabolic point operator essentially eliminates expansion shocks 
by insuring that 

0 yy = 0 (56) 

when <f> = 1. This condition forces d> to be finite. If this 
mixed difference scheme is applied to the left hand side of 
Equation (54), that is, the straining term is zero, and the term 
Uqx' ^ix' -*- s eva l uatec i with central differences then the shock 
jump condition is 

u i + U" = AxU^ x (57) 


compared to the theoretical condition 
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U 1 + U 1 = 0 (58) 
The shock jump [U^] is given by 

[U ± ] = 2U+ - AxU^ x (59) 

Equation (54) with f put equal to zero was solved using the 
mixed difference scheme described above and the pressure 
distribution around an 11% biconvex airfoil at M_ = 0.808 is 
shown in Figure 9. The base solution is a 10% biconvex airfoil 
at the same Mach number. It can be seen that the linearized 
result is smooth and that the shock is located at the wrong 
position since it is still at the base position. There are no 
oscillations in the pressure distributions in the neighborhood of 
the sonic line and consequently no apparent mechanism to 
correctly locate the shock wave. In this case the difference 
scheme, which smooths oscillatory behavior, has eliminated the 
oscillations which indicate an incorrect shock location. The 
mixed difference scheme can give plausible but incorrect 
solutions. The weakened shock jump is explained by the modified 
jump condition given by Equation (59). 

In order to obtain the correct behavior. Equation (54) is 
solved using central differencing everywhere. The result, 
without straining, is shown in Figure 10 and it may be seen that 
the solution oscillates in the neighborhood of the sonic point. 
If the straining terms are included and evaluated using Equation 
(40) then the solution is smooth in the region of the sonic lines 
can be seen in Figure 11. 

A second example, using central differencing, is shown in 
Figures 12 and 13 where the pressure distribution over a NACA0010 
airfoil is shown. Figure 12 shows the result without coordinate 
straining and Figure 13 the result with coordinate straining. 
The Mach number and base solution are the same as for the 
previous example. 
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The grid used in these computations is relatively coarse in 
the y direction (14 points) because of computer storage problems. 
This may account for the variable accuracy of the calculations. 

From the discussions above it can be inferred that the 
finite difference algorithm used to solve the TSD equation must 
be nondissipative so as to allow oscillations to occur in the 
neighborhood of the sonic line. These oscillations are removed 
by locating the shock wave in the correct position. The fact 
that the dissipative algorithm gives a plausible but incorrect 
solution should be a warning. It is imperative that the 
mathematical aspects of the differential equation being solved 
are understood before a choice of algorithm is made. 


4 . ) Linear Perturbation Equations with Constant Coefficients 
4.1 Introduction 


One of the problems that arise in the previously reported 
method for linearizing the transonic small disturbance (TSD) 
equation is that the number of shock waves in the solution is 
fixed by the base solution. Since transonic flows can have 
multiple shocks, this is a severe restriction on the method since 
the number of shocks in a solution must be known a priori. A 
second problem is the behavior at the sonic line which is 
different in a finite difference approximation of the TSD 
equation than in the mathematical formulation. In the former 
oscillations appear in the neighborhood of the sonic line, while 
in the mathematical picture an expansion shock forms. If the 
solution could be linearized about two constant velocities, one 
supersonic and one subsonic, then these problems would be 
reduced. Since the transonic solution would now consist of 
combining any number of elements of the supersonic solution into 
a subsonic solution, the shock being represented by a jump from 
the supersonic solution to the subsonic solution, any number of 
shocks could be incorporated into the solution. Also, since the 
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base solutions do not pass continuously through sonic conditions 
the problem at the sonic line would disappear. The following 
analysis suggests a method of implementing this idea. 


4.2 Basic Equations 

The integral equation for the velocity, u, can be written in 
an arbitrary coordinate system as 


u - f" - l L * + 6x o X g + 5x s I f * 


(60) 


where I T is the Prandtl-Glauert velocity, I is a field integral 
Li 

and I , I. and I, are functions from a three point coordinate 
g h r 

straining with 5x q , 5x g being straining in the x-axis and 5y g a 
straining in the y-axis . The terms I(u^/2), Ig, 1^, 1^ are 
continuous and possess continuous first derivatives. Sonic 
conditions occur when u is unity. 


Now let 


u = + u 1 

o 1 


(61) 


where u q is either 0 (a subsonic value) or 2 (a supersonic 
value) . 


Equation (60) can now be written as 


<“o + u l> (1 - V - r L + 1 ^ + Il(u o + u l> V - '("o 2 ' 2 * 


2 u 2 

+ T - * Sx s : £ + 6x o X g + S *s X h 


( 62 ) 


Let an * denote values on a sonic line, (u, + u. = 1). 



Equation (62) can be written as 


u, 


< u . 


V (i - 


u o ) = F 


2 , 

" V 2 


(63) 


where F is a continuous function. The main object is to 
determine where the solution should transfer from the subsonic 
solution to the supersonic solution and vice versa. 

. * 

Let u Q jump from zero to two at some x and from two to zero 

at x„ . 
s 

Differentiation of Equation (63) gives 


< u o + u l>x - 


V - 


U (U 

ox v o 


v ■ 


P x + u l u lx 


- u u 
o ox 


(64) 


or 


< u o + Vx <* - u o - u l> - F x < 65 > 

If the condition that u Q + u^^ passes smoothly though unity 
is imposed. Equation (66) gives the shock condition of the 
classic integral equation. 

Now linearize Equation (63) about u Q ; thus 


< u o + u i> (1 - u Q ) = F - m 2 q / 2 (66) 

where F is a linearized version of F. 

It is desired that on either side of x*, u q + u^ is unity 
(sonic). Thus, 

(u o + u x ) _ = ?* = - (F* - 2) = (u Q + u x ) + = 1 (67) 
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or 


F* = 1 (68) 

where the "+" and superscripts denote values just ahead of 

and behind x respectively. Equation (68) will give an equation 
for the shock location 

Differentiation of Equation (66) gives 


(u + 
v o 


u l>x 


(i - u o ) 


u, u = F 
1 °x 


(69) 


where a tilde denotes an average value at the point of 

* 

differentiation. Since u qx is infinite at x (it is a delta 

function) and since (u^ + u. ) is required to be bounded and F 

O -L X X 

is finite it follows that 


u* = 0 (70) 


If 



is to be finite at x 


* 


t 


* 

then since u^„ 

OX 


is infinite 


u 


* 



(71) 


Since u* is the derivative of a step of magnitude 2 it follows 

OX 

that u. is the derivative of a step of magnitude 2. Using this 

X A 

information and Equation (70) it follows that 



hence, the condition that u QX + u lx be finite leads to the result 

that the sonic line occurs at x . In other words the 

* 

linearization of the equation makes the sonic line at x . Since 
x can be put anywhere it follows that the sonic line location is 
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arbitrary. It is only accurate if the (neglected) nonlinear 
terms in Equation (64) cancel. That is if 


~ * 



~ * 



I* x (u^/2) 


(73) 


From Equations (62) and (63) 

V 1 - V - J L + "W + t 1 *^) - % + u o /2 i 

+ 6x s 1 i + 5x o X g + S *s 

= F (74) 


where F is continuous. Expanding Equation (74) as a Taylor 
series about x gives 


| /V A 

u 1 = F - AxF x 

A A 

U~ = -(F + Ax F x ) (75) 


Thus 






A ^ A ^ 

(-2AxF x ) (-2 F ) 
2 2Ax 


F x F = F x (76) 


* 

where it has been noted that F = 1 . 

It can be seen that u ^ u i x actually a first order 

quantity. Hence, a linearized version of Equation (73) is 

F* = 0 (77) 

Consider now Equation (66) at values just ahead of and behind the 
"sonic" line and denoted by the superscripts "+" and 
respectively. On addition. 


o n 


( 78 ) 


\ [<U Q + u 1 ) + + (U Q + n x ) ] = \ [F + + F ] - 1 
or, since F is continuous. 


F* = 0 (79) 

Finally, let the shock and the sonic line coincide at the point 

(x ,y ) (the "top" of the "shock" in u ). This gives 
s s o 

< e *o + Vy = < 5 x s + V'y (80> 

s s 

Thus the governing equations are Equation (66) with Equation (79) 

applied at x*. Equation (69), Equation (77), and Equation (80). 

These are sufficient for solving u, at N discrete points, 5x q , 

frx on each grid line and 5y . 
s s 

As a first approximation the Sy„ term is neglected and only 

s 

Equations (66), (69), (77), and (79) need be solved. 

4.3 Possibility of Deriving a Closed Form Solution 

In the classic strained coordinate theory 1^, 1^, 1^ are 
evaluated using the base velocity, u q . However, since u q is 
piecewise constant it is instructive to use the unknown (u q + u^) 
in this evaluation. 

Consider the equation 

(1 - u ) <b + <b =0 (81) 

v o ' r xx Y yy v 


Since u q is constant. Equation (81) can be written as 
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^x'x' + ^yy ^ u o^x'^ + ^ u o^^x^x'^x' 


f x < X - u o)t< 1 +f x)^x'lx' 


where 


x' = x + f(x) 


In the general coordinate straining, x' is given by 


x = x' + x 1 (x' ) 


*x - < x + £ x>*x- 


It follows from Equations (83) and (84) that 


* + ^x 1+x, 


and hence Equation (82) becomes 


^x'x' + ^yy ( u o^x')x' + ^ X 1 , ^ u o^x^ , 


^1+X- ^ ^ V ^xx' 

x' 


If the terms involving X 1 , are linearized, then Equation (87) 


becomes 
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*x'x- + *yy - <Vx-> , + t x l x ' (1 “- V *x> , 


+ X 1 , <1 - V*XX' 

X 


( 88 ) 


Equation (88) can be regarded as a piecewise application of the 
subsonic and supersonic Prandtl Glauert equations in a strained 
coordinate system. Thus in the subsonic domain, Dsub, the 
equation is 

^x'x' + ^yy = t x lx' ^x'^x' + x lx' ^x'x' (89) 


and in the supersonic domain, 
~^x'x' + ^yy = “ [x lx' 


Dsup, 



the equation is 
X 1 x'^x'x' 


(90) 


where 6 is some suitable (zeroth order) approximation to <p , . 
These equations could be solved piecewise analytically in 
specified subdomains, Dsub, Dsup, with the matching condition of 
Equation (72) together with the condition 

v + = v“ (91) 

This second condition is implicit in the derivation of the 
integral equation. Equation (62). Both of these conditions are 
sufficient to choose a strained coordinate system but, as in the 
case of the integral equations, the sonic line will always be at 
x . It is not clear how a condition similar to Equation (73) can 
be implemented the differential form of the equation, and because 
of this a purely analytic solution of Equations (89) and (90) is 
probably useless. 


II 

II 

II 

I 

II 

II 
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I 
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4.4 Results 


Some results were found using the integral equation theory 
with a triangular hyperbolic region. The height of the region 
was kept fixed, that is , no y straining is used. The hyperbolic 
domain closely resembles that for the 10% biconvex airfoil at 
M = 0.808. The pressure distribution over a NACA0010 at 
M_ = 0.808 is calculated and shown in Figure 14. The straining 
terms are evaluated using the base solution. A second result is 
shown in Figure 15; in this case the straining terms are 
evaluated using an estimate of the final, smooth solution. The 
results are very sensitive to small changes in I L , and it is not 
suggested that this is a viable computation procedure. However, 
considering the size of the perturbation from the base solution, 
the accuracy of the results is surprising. 


5. Perturbation of the Euler Equations 
5 . 1 Introduction 


In section (3) a method of locating the shock wave in the 
linearized TSD equation that is suitable for finite difference 
solutions is derived. If these ideas are extended to solutions 
of the linearized Euler equations in this section, it is found 
that the shock wave can be located correctly but only if the odd 
and even points are separated in the solution. 

5.2 Basic Equations 

The Euler equations can be written in general coordinates 

(£> V) as 


A A 



( 92 ) 
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where 



>U 


pV 

-1 

pu u + 4 x p 

* -1 

puV + rj x P 

E = J 

pv U + £ y p 

; F = J 

pW + »/ y p 


-(pe + p)U. 


.{pe + p)V. 


(93) 


where (x, y) is a cartesian coordinate system with velocity 
components u, v in the x and y directions respectively, p is 
pressure, p is density and e is given by 


12 

e = e i + 2 q 


(94) 


where e ^ is the internal energy and g is the total velocity. U 
and V are given by 


u = £ u + e v 
s x s y 

V = » 7 x u + rj y v 

p, u, v are normalized with respect to freestream values, denoted 

2 

by the subscript oo, and p and pe are normalized by /) tt q^. J is 
the Jacobian defined by 


(95) 


J = - Vx 


(96) 


Equation (92) can be written as 


a aaa a aaa 

[A (Q) Q] + [B (Q) Q] = 0 

i 

i where 


(97) 


Q 



[P, 


n T 

you, pv, e] 


(98) 


I 
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and 


> > 
II 

0 

~ uU + £ x 0 2 

U + (2 - 7 )C x u 

Cy U - (7 - !)C X V 

0 

(7 “ 1)C X P 


-vu + £/ 2 

€ x v - (7 - l)£ y u 

U + (2 - 7KyV 

(7 " l)C y P 


- aj 

C x a i - (7 " !)uU 

e y a i - (7 - 1)VU 

7 U 





( 100 ) 

A 

B = 

0 

-uV + V 2 

v + (2 - 7)7 x u 

V Y 

V y u - (7 - l) 7 x v 

0 

(7 - !)7 X P 


-vV + r}y<p 2 

V - (7 - l)7 y U 

V + (2 - 7 ) 7 y v 

(7 ~ l)»7yP 


V[(f> - a x ] 

7 x a i - < - l)uv 

7 y a i - (7 - 1 )W 

7 V 

where 





a l 

= 7 e//> - 0 2 ; 0 2 = (7 

- 1 ) (u 2 + v 2 )/2 

( 101 ) 

and 


£ = 4 £ = 

9Q 

A 

OF 

A 

9Q 

( 102 ) 


For a nonlifting airfoil, the numerical solution of Equation 
(97) requires the following boundary conditions. 

On the far field boundary, freestream values are imposed on 
all variables; this is valid strictly only for isentropic flows. 
On the body surface, rj = 0, 
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V 


9£ 

drj 

9U 
9 V 

9e 

dr) 


0 

0 

0 


( 103 ) 


These boundary conditions impose zero vorticity on the airfoil 
surface, again valid only for isentropic flow. Although the 
imposed boundary conditions leave something to be desired, it is 
expected that useful insight into the solution of the linearized 
Euler equations can be gained. 


5.3 Linearized Euler Equations 


The linearization of Euler equations used in the present 
analysis is designed to be similar to that used for the TSD 
equation and, in fact, reduces to the perturbed TSD equation if 
the correct limiting processes are used. 

Equation (97) is written in a general coordinate system so 
if the x coordinate is strained such that 


x = x' + ex 1 (x' ) 


(104) 


where (x',y) is the strained coordinate system with a straining 

parameter, 6 , then the effect can be felt only through the 
xs 

Jacobian, J. 

The coordinates are strained as follows 


£(x, Y) = £( x ', y) 
»7(x, y) = *7( x , y) 


( 105 ) 
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such that 


~ ^x' ex lx' ^x' 

^x ~ V + €X lVx' 


( 106 ) 


The relations of Equations (105) and (106) are similar to those 
used in the TSD theory. The general variable Q is expanded as 

Q = Q q + eQ 1 (107) 


where 


Q = [p, pu, pv, e] T 


(108) 


A A 

The Jacobians A, B can be written as 


A = A o (Q o ) + 6A l( Q 0 ) 
B = B o (Q o ) + 6B l( Q 0 ) 


(109) 


where eA 1 and eB 1 are the perturbations of A and B respectively 
due to a change in £ , rj only. 

A A 

Finally, the Jacobian J is expanded as 


J = J + eJ. (110) 

O 1 

Substitution of Equations (107), (109), and (110) yield the 

following perturbation equations 
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(Ill) 


5? [A 0 (Q 0 )e 0 ] + 8^ [B 0 (Q 0 )Q 0 1 - 0 


ft A Q A ^ Q A A a W 1 A A 

a| tA o^l ] + drj [B o Q l^ = 9f [ j“~ A o Q o ] + 8*7 B o Q o ] 



A A a A A 

l A l Q ol - drj ("W 


( 112 ) 


It should be noted that 

x lx' ^x' *o 
x lVx' § o 



(113) 


where A and B are cartesian Jacobians and 
o o 



(114) 


The boundary conditions for the zeroth approximation. Equation 
(111), are the same as those for Equation (92). 

In the particular formulation of the Euler equations used 
here a change in freestream Mach number, M w , is felt only through 
the far field boundary condition on p and e. For example. 


P 


00 



«. = \ + 1/11(7 - 1>M»] 


(115) 


In order to simplify the algebra a simple test case is chosen. 
The test case is a flow with a perturbation in freestream Mach 
number, the geometry remains fixed; this allows the same grid to 
be used for both the zeroth and first approximations. A second 
simplification is that the airfoil is symmetric and at zero angle 



of attack; this allows the use ' of symmetrical boundary 
conditions, thus reducing the computer requirements. 

The boundary conditions for the perturbed, primitive 
variables are therefore as follows . 


At the far field boundary 


(P u >! 

(P v >! 

e l 

P 1 


0 

0 


7(7 - 1)M, 
1*1 


(116) 


where e is a perturbation parameter. The perturbed flow is one 
where the freestream Mach number, M^, is given by 


M l = **„/<! + e> 


(117) 


On the airfoil surface, denoted by rj = 0, boundary conditions 
similar to those in Equation (103) are used. 

On the symmetry planes the boundary conditions 


9^1 


av 


i 




3e 


1 




0 

0 

0 

0 


( 118 ) 


are used. 


For the base solution. Equation (111) is solved using 
central differencing with added dissipation. This is to ensure 
that the base solution, Q q , is smooth except at shock waves. The 
linear perturbation equation. Equation (112), is solved using 
central differences with no added dissipation. At the boundaries 
first order derivatives are approximated by first order 
differences. Although done mainly for computational use, this 
has the effect of weakly coupling the odd and even points in the 
solution. Equation (112) can be written in the discretized form 

MQp = B + S (119) 

where L( ) denotes the difference operator, approximating the 
left hand side of Equation (112), B denotes the boundary terms 
applicable to Equation (112) and S denotes the terms due to the 
straining function. A simplified form of Equation (119) is 
obtained by neglecting the straining term, S; thus 


L(Q X ) = B (120) 

The solutions of Equation (119) or (120) are obtained by direct 
inversion of the operator L. Thus Equations (119) and (120) 
become 

Q ± = L -1 (B) + L -1 (S) (121) 


and 


Q 


1 



(B) 


( 122 ) 


respectively. The computations in this analysis are performed on 
a grid of 40x12 for the half plane. This is the maximum grid 
size for direct inversions. 
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The pressure coefficient found from Equation (111) is shown 
in Figure 10, and it can be seen that it is typical of those seen 
for transonic flow. The density found from Equation (120) is 
shown in Figure 17, and it can be seen that it oscillates 
considerably. This oscillation is typical of central difference 
approximations when no dissipation is added. 

In general a central difference approximation will decouple 
the odd and even points, and, although the treatment of the 
boundary conditions allows some coupling between odd and even 
points, it is instructive to look at a solution composed only of 
odd or even points. In Figures 18a and 18b the density variation 
for the odd and even points respectively are shown. It can be 
seen that the solution in the neighborhood of the sonic line is 
smooth except for a "spike;" the oscillations near the shock 
wave are due to the presence of a "nonphysical" sonic line in the 
shock capture region. Hence, it can be surmised that the same 
phenomena apparent in the TSD results exists for the "decoupled" 
linear Euler equations, that is, the solution is smooth except 
for an oscillatory behavior in the sonic region. In analogy with 
the TSD result the new shock location is enforced by requiring 
that the decoupled solutions are smooth at the sonic line. This 
is accomplished in the same manner as for the TSD equation, 
namely by enforcing the condition that 



(123) 


at the sonic line. Thus, the solution procedure is to solve 
Equation (121) and determine the parameter 6x in the term L -1 (S) 
by enforcing Equation (123). The complete solution for Q is 
determined by 


Q(x, y) = Q 0 (x', y) + eQ 1 (x', y) 
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5.4 Results 


The pressure distribution around a NACA0010 airfoil with 
M ro = 0.82 is shown in Figure 19 with a base solution of the same 
airfoil at = 0.8. In Figure 19a the odd points are shown, and 
in Figure 19b the even points are shown. It is not clear why the 
even points should give a more accurate result than the odd 
points . 

In Figure 20 the pressure distribution around a 10% biconvex 
airfoil at = 0.828 is shown with a base solution of the same 
airfoil at = 0.808. The odd and even results are shown in 
Figures 20a and 20b respectively and again it may be observed 
that the even points give a more accurate result. 

There is a general problem with these Euler computations 
which is the relatively low number of grid points . The maximum 
number of grid points is determined by the maximum in core 
storage of the CRAY X MP machine to allow a direct inversion of 
the operator L( ) . It is expected that the grid error in the 
decoupled points is considerable. 

5.5 Concluding Remarks 

A linearized form of the Euler equations has been developed 
and a method of locating shock wave derived. The key to the 
method is the decoupling of odd and even points and the removal 
of the characteristic oscillations at the sonic line. 


6.) Analysis of the Modeling of Separated Flow Using the Euler 
Equations 


6.1 Introduction 

In recent years solutions of the Euler equations have been 
used to model the flow around aircraft or missiles. The main 
advantage of the Euler equations over simpler sets of equations 
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or equations such as the potential equation is that the Euler 
equations can model nonisentropic, and hence rotational flows. 
In certain cases, especially missile flows, the ability to 
capture contact discontinuities and hence model separated flow 
has been an important feature of the Euler equations. Such flows 
have several sources of vorticity. Vorticity is generated from 
curved shock waves and in the viscous boundary layer on bodies . 
If the boundary layer flow separates from the body, a large 
vortical region is produced which can be modeled by the Euler 
equations if the necessary vorticity can be produced. For the 
Euler equations, vorticity can be produced by curved shocks, by 
certain boundary conditions, or by numerical discretization 
error. The numerical error is a somewhat nebulous quantity which 
seems to be controllable only by careful smoothing and 
optimization of the computational grid. Boundary conditions for 
modeling separation usually involve specification of the 
velocities at or around a specified separation point. If the 
separation is from missile fins then the separation point is 
taken to be the fin edges; if body separation is to be modeled 
then some empirical separation criterion is used to determine the 
separation point. Also, because of the vorticity generated by 
curved shocks and numerical error, the Euler equations can have 
solutions which contain separated flow regions even without the 
specification of a separation point. 

In this section various aspects of modeling separated flows 
on smooth surfaces using the Euler equations are examined. It is 
found that "Euler separation" can be radically different in 
behavior than real, physical separation. It is found also that a 
separation point cannot be specified arbitrarily on the body 
without further modification to the Euler equations. It is 
suggested that numerical treatments of the boundary condition may 
not be consistent with the Euler equations. 

6.2 Analysis 

The Euler equations are examined to determine the conditions 
under which separation is possible in the sense that the flow 
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tangential to a surface can have a stagnation point 
vorticity equation. 


U - V = w 

y x 


can be written in terms of a stream function ^ such that 


*xx + I'yy • p u + V - W 


where u is the vorticity and 


*y’PV- PM * X * -P V 


The density is given in terms of the entropy, S, 
velocity components U, V by 


P = Poo { 1 + il 2 1) C 1 - U 2 - V 2 ]} 7 1 exp(-S/R) 


The vorticity is given by Crocco's theorem as 


t as p # as 

q 3n R^q 3n 


where R is the gas constant 


q 


(U 2 + V 2 ) 1/2 


and 


a_ 

3n 


The 


( 124 ) 


( 125 ) 
and the 


( 126 ) 


( 127 ) 
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is the derivative in the direction normal to the streamlines. 
Note that 


L = U L. V L 

an u ay v ax 


Along a streamline the entropy is constant and hence 




The pressure is given by 


P = P w {1 + J ^ kl M 2 [1 - U 2 - V 2 ] } 7 1 exp( -S/R) 


Using Equations (126), (127), (128), and (129), Equation 

can be written as 


- *vv = p« Q 7/7 " 1 exp(-s/R) <B || - V If) 


xx r yy 


Rq 


3y 3x ; 
1 


+ P* IQ 7 1 «xp(-S/R)] U - p w [Q 7 ixp(-S/R)]V 


where Q is a function of U and V and is given by 


Q = [1 + M 2 (1 - U 2 - V 2 


)} 


Equation (131) can be written as 


^xx + *yy = f < x ' y> 


where f(x,y) is the right hand side of Equation (131). 


(128) 

(129) 

(130) 
(124) 

(131) 

(132) 

(133) 
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The problem under consideration is the flow around a thin 
body, given by 


y = 

■*su< x > 

rv 

0 



(134) 


7si (x) 

VI 

>i 

0 

where the subscripts u and 1 

denote the 

upper and lower surfaces 


of the body. The body may be at an angle of attack, a 

To some degree of approximation it can be assumed that the 
body boundary condition can be applied on a slit y=±0. The first 
approximation is thin airfoil theory which gives 



r > 

*su< X > 

K u oojy =!t0 

7sl< x >. 


(135) 


Equation (133) can be written in integral form using Green's 
identity to give 


J D /K(x,£;y,y)f(£,T/)dD = - f {K(x,f ;y,? 7 )|^ (€, 1 J) 

J “ / ^ n in in i n 


c.+c +c +c 

1 w s 00 


- tU.D )dc 


(136) 


where C is the boundary of some domain D and v is the unit inward 
drawn normal to C. The domain D excludes all regions of the flow 
where second derivatives of ^ or K do not exist, such as the body 
(represented by y = ±0), wakes, shock waves, and contact 
surfaces. In Equation (136) is that part of the boundary 
surrounding the point (x,y), C w represents the boundary 
surrounding the x axis, C the boundary surrounding the shock 
waves or contact surfaces and is the far field boundary. A 
sketch of the boundaries C- , C , C , C are shown in Figure 2; to 
avoid confusion with the entropy, the domain S in Figure 2 is 


i 
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denoted by D in the present analysis. The kernel function in 
Equation (136) is given by 


1 

K(x,f;Y/»7) = 2?r ln {< x ~€> 2 + (Y~V) 2 } 2 


(137) 


For subsonic flow the Euler equations are elliptic and hence 
boundary conditions must be specified on the entire boundary. 
For transonic flow the flow outside of a finite domain is 
subsonic so that the same conditions must be satisfied on the 
outer boundary, C^. Equation (136) is similar to the integral 
equation for potential flow, and it may be deduced (Ref. 6) that 
the line integral over will contribute a finite velocity on C M 
as R-»-<» if 


if> = A InR + B 


(138) 


where A and B may be functions of the angular coordinate 6 . 

If this relation is used then the line integral over 
becomes constant . 


On performing the various limiting operations as r-»-0 and R-*-°o 
Equation (136) becomes 


*( x ,y) = l < K o A o *n - %V }d S 

COS(l/,»7 ) 

- /c-fVV- + K s cos („,£„ ) [ty_> d *7 

s s 

. 9K 9K cos(i /,« ) 

* /c co.,,./, > d " 


+ J D /Kfdfd»j + constant 


( 139 ) 


where for a function g(£,»7) 
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(140) 


A 0 g = g(£, + 0) - g(C, - 0) 


and 


K q = K(x,£; Y/O) 
K s = K < x '£ s ; Y ' V) 


(141) 


£ and rj are the coordinates of any surface of discontinuity and 
s s 

cos (i / ,ti ), cos(i/,f ) are the direction cosines of the normal to 
s s + 

these surfaces. [ ]_ denotes a jump across the surface. 

If the jumps in derivatives of ip , such as and ip , are 

x y 

finite then rp must be continuous since a jump in ^ would give an 
infinite derivative. Hence 

= 0 (142) 

The integrand of the first integral over C in Equation (139) 

b 

becomes 

COS(l/,t) ) - 

V'V- - IV- COS(l/,{ s ) } = 'Mtl- OOS <‘''«s> (143) 

fs* 

where q t is the velocity tangential to the curve C g . Hence 
Equation (139) can be written as 

1 *<x,y> * /” 1K oV> , - K oy A o' >,d$ 

- I c K s [/9q t ]%ec(i/^ s )dr/ 
s 

+ / /K fdfd rj + constant (144) 



Equation (144) can be differentiated with respect to y to give 


t b = pU-p U = r °°{ K A it - K A it} df 
' y r ^00 00 J O x Oy O r Yj OTjy O r J s 

" f c K sv [ P q t ] - sec(I/ '€s )<ir l 

+ ; D /K y fd«d, 

By integration by parts it can be shown that 

/o K o W V d f - /o K oxV 5 d « 

loIKytWV - /“ K oA Fd « - 


where 

F = I v f(Z,ri')av' 

Thus Equation (145) can be written as 

pv-PA - - /Xx A oV« + 'oW, - A o F]d « 

“ IC K S y lP ^t]- Sec ( 1 ''« s ) d ’7 
s 

- 

On differentiation of Equation (144) with respect to 
performing a similar integration by parts 


(145) 

(146) 

(147) 

(148) 

(149) 
x and 
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( 150 ) 


pV = - / K {A - A F}df + / K A df 
r J o ox 1 o r iy o J ^ J o oy o r ry * 

- Ic K sx (Mt]- SeC<l/ 'fs )d ’’ 

s 

- V K x, Fd ^ 


If the limit as y ■*■ • o is taken of Equations (149) and (150) the 
following equations result 


y = ± 0 1 T - ^ U »“ - C 

4 X f 

- Lim. { J -fK [F(£ ,rj) - F(x,±o)]d£d »7 
y -*±0 * ' 



sy 


[pq t ]* sec < l/ /C s ) d *7} 


(151) 


- {(pV) 


y = 


+ A o<' v > 

,00 

> = /« 

[A o <p D >-A o : 

±o 2 

* o 

x - e 


+ Lim. 



y-»±o 

c s 



- V K x, 


-d£ 


(152) 


It can be seen then that on y = ±o Equation (149) gives 
only the symmetric part of pU; the asymmetric part is given by 
Equation (152) in terms of the symmetric part of (pV)y _ Q . 

It should be noted that in both Equation (149) and Equation 

(150) the integral over C may contribute a discontinuity 

s 

depending on the orientation of C to the axis. For example, if 

s 

C is in the direction of the v axis (a normal shock wave 
s 

perhaps) the integral over C in Equation (149) is continuous. 

s 

Since the other integrals are continuous (Ref. 6), pU is 
continuous through a normal shock; this is the correct Rankine- 
Hugoniot solution. However, in Equation (150) the integral over 
C provides the term 



- 2 [P V l^ s 9 n ( x s -x ) 

where x is the location of the shock. This term ensures that 
s 

Equation (150) is consistent on both sides of the shock. The 
actual value of [pV] + can be found from the jump relations, 
either for a shock wave or a contact discontinuity. 

From Equations (126) and (131) it can be seen that p and f 
(and consequently F) are functions only of U and V and the 
entropy S. The entropy is given by Equation (129) and the 
velocities just ahead of the shock wave. Thus S is a function 
only of U, V. Hence, if (pV) y = ±Q is defined as a boundary 
condition Equations (149) and (150) are equations for U and V 
throughout the flow field. Because of the complexity of the 
relations for p , S and F it is not easy to determine the number of 
possible solutions. However, it is safe to say that if pV is 
specified on the body then pU and pV are determined completely in 
the flow field; pU is also determined on the body. In 
otherwords, U and V cannot both be specified on the boundary; 
this statement requires further discussion. 

It has been the practice in some problems to model separated 
flow using the Euler equations. In these methods an empirically 
determined separation line is specified. The separation line is 
such that the flow on at least one side of the resulting contact 
surface on the body surface is a stagnation point. In effect the 
method is specifying both U and V at this point and from the 
above argument this overspecifies the boundary conditions for the 
Euler equations. The necessary additional degree of freedom for 
these techniques could be introduced by an additional entropy or 
vorticity source which would mimic the absent boundary layer 
generated vorticity. 

A second point that arises is the use of extrapolation of 
the tangential velocity to get a value on the surface, for 
example, by the use of boundary conditions of the type in 
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Equation (103). These conditions specify the tangential velocity 
on the surface in contradiction to the above argument. 
Consequently, solutions employing this device are not solutions 
of the Euler equations, although the error is difficult to 
estimate, and the approximation may be adequate in an engineering 
sense. 


The asymmetric or “lifting" part of pU is given by Equation 
(152). This equation is similar to one in classic aerodynamic 
theory and has at least one eigensolution such that Equation 
(152) can be written as 


+ A o<^ V > > = r t A o<' U >- A o F - AG <OJ 

y=±o 2 ° » - * a $ 


x - i 

+ Lim. {-jK sx [pq t ]^sec(i/,f s )d»? 


y+±o c c 


- W d£d,} 


(153) 


where 


G(0 


1 

ri/te <i - o i 2 ; 

to ; Oi 


(154) 


and A is an arbitrary constant. An additional condition is 
necessary to determine the constant. A general solution of 
Equation (153) is 


A o (pU) - A q F = AG ( x) + H(x) (155) 

where H(x) is the "basic" solution. Since any choice of A other 
than zero would make (A Q (pU) - A Q F] infinite at the trailing edge 
the logical choice is 
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A = 0 


(156) 


It is important to note that this condition, the classic Kutta 
condition, must be imposed either explictly or implicitly in any 
solution procedure since, in general, the Euler equations can 
give an arbitrary lift. 

As a final point it can be stated that if a separation point 
in the Euler equations is defined as a stagnation point, then for 
a symmetric body Equation (151) can conceivably give a solution 
that will make U zero at some point on the surface. The actual 
point is determined by the body geometry and by the entropy 
generated in the flow. For a lifting flow a similar situation 
arises but some measure of control on the location of the 
separation points is given by the arbitrary constant, A, 
discussed above. 

6.3 Additional Comments on "Euler" Separation 

One of the jump conditions across a vortex sheet is that 

Ap = 0 (157) 


Equations (130) and (157) then give 
{1 + 13 z. 1 ) [1 - q +2 ]} exp(-S^Cp) 

= U + ^ 2 ^ M » t 1 “ q" 2 ]} exp(-S"/c p ) (158) 

where "+" and signs denote values on the upstream and 

downstream sides of the separation point respectively. 

If the velocity on the downstream side is zero, then 
Equation (158) gives 
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TZ 


{[1 + 


m;] [exp(AS/c p ) -1]> (159) 


q 


£ 

( 7 - 1 ) 


vi- u. 

2 


where 


AS = S + - 


S 


(160) 


For a real value of q+ to exist 


exp (AS/c p ) S 1 


(161) 


For AS equal to zero, that is isentropic flow, 

q + = 0 (162) 


In other words, if separation occurs in isentropic flow, the flow 
has a stagnation point at both sides of the separation. 

From Equation (161) 


S + < S (163) 

If a shock wave is upstream of the separation and produces 
entropy then Equation (163) cannot be satisfied. Re-examining 
Equation (158), it is found that this implies that at the 
separation point the stagnation point should be on the upstream 
side of the separation point, in other words the separation 
streamline would be inclined in the upstream direction. This 
flow feature does not appear in experiment thus indicating some 
of the peculiarities of separation phenomena in the Euler 
equations . 



Since the flow is attached before separation it follows that 
for a convex surface 


9n 


> 0 


(164) 


where n is the direction normal to the surface. 

Using Equation (130), Equation (164) becomes 

-1 M »3 " (1 + 1SL T 1 M l (l-q 2 )} (S/R) * 0 (165) 


or 


0g 

9n 


< - 


<— r> f 1 + l 2 r 1 m « ( i -< j 2 > j Ik < s / r > 

7 M «q 


(166) 


If the flow is isentropic then 9q/0n is negative in contradiction 
to 9q/0n for a viscous flow which is positive. However, if 
9(S/R)/0n is negative, as it will be for a shock generated 
entropy, 9q/9n will be positive and have the same sign as for a 
viscous flow, thus giving some realism to the flow model. 

Separation is characterized by the vorticity in the boundary 
layer . For a turbulent boundary layer an approximation to the 
vorticity is given by 


w, - 
b 


93 ~ a 

9n " g 


0_ 
e 9n 


(£> ? - 


15 


(b 


-6/7 


(167) 


where w b is the vorticity due to the boundary layer, 5 is the 
boundary layer thickness, and q e is the velocity external to the 
boundary layer. From Equation (167) 


Ui 


b 



(168) 
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The vorticity generated by the shock wave is given by Crocco's 
theorem. Equation (129), and can be written as 


1 af axs/Rl 

E 7 q 9n 


(169) 


where u & is the vorticity generated by the shock wave and a is 
the speed of sound. Hence, 



7 6 

2 

1 % 


3(S/R ) = 
9n 


7 5 

gM 2 


k < S/R > 


(170) 


where M is the local Mach number for the Euler calculation and q 

has been identified with q . 

e 

The entropy rise through a weak shock is given approximately 
by (Ref. 7) 


0(S/R) 

9n shock 




( 7 + 1) 2 



- l) 2 ^ 


0 M 1 

dn 


(171) 


where M^ is the Mach number just ahead of the shock. If it is 
assumed that 


k < S/R > 


c 9_ ( S/R) 

3 n shock 


(172) 


then Equation (170) gives 


U) _ 

_E < 7 C 5 

u , " .,2 

b 7 M 


— 4 7 

( 7 +l ) 2 



- i) 2 m x 


9M x 

3n 


shock 


F 


(173) 


This equation indicates that the vorticity generated by the shock 
is much less than that generated by the boundary layer for weak 
shocks but that as M^ increases, the inviscidly generated 
vorticity can dominate. 



A crude indication of the magnitude of is given as 

follows: Let C be unity, M - 1, 5 =0.01, (the boundary layer 

thickness). Finally, assume the shock length is unity and that 
the variation of M^ along the shock is linear. Hence, 

0M 

g^T “ -< M 1 " X > < 174 ) 

The variation of the function F with M^ is shown in Table 1. 


M 1 

F 

1.0 

0 

1.2 

0.0023 

1.4 

0.0251 

1.6 

0.1136 

1.8 

0.3512 

2.0 

0.8750 

2.2 

1.8924 

2.4 

3.7007 

2.6 

6.7092 


Table 1 

It can be seen that the overall contribution of the shock 
generated vorticity increases rapidly with the strength of the 
shock, indicating that at the higher supersonic Mach numbers the 
vorticity field is dominated by the shock induced vorticity. 


7 . ) Concluding Remarks 

This work is concerned with developing methods for solving 
transonic flow problems using linearized forms of the TSD 
equation and the Euler equations. Methods have been developed 
using both integral equation methods and finite difference 
methods. A key element in both the TSD and Euler models is the 
use of a strained coordinate system in which the shock remains 
fixed. Additional criteria are then developed to determine the 
free parameters in the coordinate straining; these free 
parameters are functions of the shock location. 
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From an integral equation analysis it is found that the 
shock wave is located by ensuring that no expansion shocks exist 
in the solution. If a numerical solution is obtained, this 
expansion shock is represented by oscillations in the solution 
near the sonic line; the correct shock location is determined by 
removing these oscillations. In other words, the oscillations 
are a crucial element in locating the shock wave. This 
conclusion is true for finite difference solutions as well as 
integral equation solutions. Finite difference algorithms 
frequently have dissipation to remove oscillations caused by 
numerical error at shock waves. It is shown in the present work 
that such algorithms can remove the oscillations at the sonic 
line, thus eliminating the mechanism for locating the shock wave. 
An example is given of a plausible result which has the shock in 
the wrong location. The important point is that algorithms 
should reflect the mathematics of the equations. Since in many 
cases the mathematics are unknown, it is possible that commonly 
used algorithms could lead to nonphysical results. 

A second, major, study is that into the ability of the Euler 
equations to model separated flow. The investigation shows that 
the correct boundary condition is that velocity normal to the 
solid body should be zero; all other flow variables can be 
obtained from the resulting equation. Thus it is not consistent 
to specify a separation point on the body since this imposes an 
inconsistent value of the tangential velocity. The one exception 
to this statement is for a lifting case where a Kutta condition, 
which imposes a separation point, is necessary to close the 
solution. As a final point, it is shown that "Euler separation" 
does have some nonphysical features. Again, it must be stressed 
that more study of the mathematical nature of the Euler equations 
is necessary to prevent the use of algorithms that are 
inconsistent with the differential equations. 


- 58 - 



REFERENCES 


1. Jameson, A. and Baker T. J. : Solution of the Euler 

Equations for Complex Configurations. AIAA Paper 83- 
1929, 1983. 

2. Nixon, D. and Liu, Y. : Concepts for Radically 

Increasing the Numerical Convergence Rate of the 
Euler Equations. Nielsen Engineering & Research, 
Inc. Report 326, 1984. 

3. Nixon, D. : Perturbation of a Discontinuous Transonic 

Flow. AIAA Journal, Vol. 16, No. 1, 1978. 

4. Nixon, D. : Calculation of Unsteady Transonic Flow 

Using an Integral Equation Method. AIAA Journal, 
Vol. 16, No. 9, 1978. 

5. Stahara, S. S.: Operational Manual for Two- 

Dimensional Transonic Code TSFOIL. NASA CR 3064, 
1978. 

6. Nixon D. : Calculation of Transonic Flows Using 

Integral Equation Methods. Ph.D. Thesis, University 
of London, 1976. 

7. Liepmann, H. and Roshko, A.: Elements of Gas 

Dynamics, Wiley, 1957. 


- 59 - 




6 





















Figure 6 
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Figure 7.- Pressure Distribution Around a NACA 64A006 

Airfoil, ML = 0.875 
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Figure 8.- Pressure Distribution Around a NACA O' 
Airfoil; M = 0.808 with Double Straining 
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Figure 9.- Finite Difference Solution with Mixed 
Differences . 





Figure 10.- Finite Difference Solution with Central 
Differences (11% biconvex airfoil). 






















DIRECT SOLUTION 
PERTURBED SOLUTION 



Figure 16.- Base Solution for Euler Equations (NACA 0010 
Airfoil ) . 







DIRECT SOLUTION 
PERTURBED SOLUTION 



(a) Odd Points 

Figure 18.- Solution of Linear Euler Equations by Central 
Differences with no straining. 


77 







DIRECT SOLUTION 


PERTURBEO SOLUTION 
(no straining) 



(a) Odd Points 

Figure 19.- Solution of Linear Euler Equations 

Coordinate Straining (NACA 0010 Airfoil) 


with 
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DIRECT SOLUTION 

PERTURBED SOLUTION 
(no straining) 
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(b) Even Points 


Figure 19.- Concluded. 
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DIRECT SOLUTION 

PERTURBED SOLUTION 
(no straining) 



(a) Odd Points 


Figure 20 


Solution for an 11% Biconvex Airfoil. 






